;+
; NAME:
;       CHECK_GRS_HI
;
; PURPOSE:
;       Check the KDA assignments of BGPS objects associated with GRS
;       clouds by way of HISA maps.
;
; CATEGORY:
;       EMAF
;
; CALLING SEQUENCE:
;       CHECK_GRS_HI [,ANGSEP=angsep]
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       ANGSEP -- [arcmin] Angular separation between BGPS peak flux
;                 and GRS cloud center to be used for assigning GRS
;                 association [Default = 10'].
;
; KEYWORD PARAMETERS:
;       NOPLOT -- DO not create the .eps plots for each source, merely
;                 run the remainder and make plots at the end.
;       UPDATE -- Force re-run of ADD_GRS_INFO.
;
; OUTPUTS:
;       NONE  [Series of EPS files, one for each source]
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  07/12/12, TPEB -- Initial version.
;       Modified: 10/30/12, TPEB -- Codified, cleaned up, and moved to
;                                   emaf/ directory.
;       Modified: 11/01/12, TPEB -- Added /NOPLOT and /UPDATE
;                                   keywords.
;       Modified: 12/03/12, TPEB -- Changed default ANGSEP to 10'
;
;-

PRO CHECK_GRS_HI, ANGSEP=angsep, NOPLOT=noplot, UPDATE=update
  
  ;; Parse
  IF n_elements(angsep) EQ 0 THEN angsep = 10.
  plot = ~ KEYWORD_SET(noplot)
  update = KEYWORD_SET(update)
  
  ;; Load the EMAF structure
  restore,'./emaf_paper/bgps_emaf.sav',/ver
  mfn = './emaf_paper/save_files/morph_distances.sav'
  IF ~ FILE_TEST(mfn) THEN run_morph_matching,/drivemorph,/silent
  restore,mfn,/ver
  
  ;; Check that the GRS elements of EMAF have been filled in...
  IF emaf[0].grs_angsep NE angsep OR update THEN BEGIN
     add_grs_info, ANGSEP=angsep
     undefine,emaf
     restore,'./emaf_paper/bgps_emaf.sav',/ver
     n = n_elements(emaf)
  ENDIF  
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  ;; Look at only the subset of sources with GRS association
  ind  = WHERE(emaf.grs,nind)
  cnum = emaf[ind].cnum
  emaf = emaf[ind]
  dml  = dml[ind,*]
  dbar = dbar[ind,*]
  print,nind
  
  ;; Declare directory shortcuts
  root  = './local/output/grs_hi'
  post  = './local/postage/'
  gpost = './local/glimpse_postage/smooth33arc'
  
  ;; Also, write out text file for handwriting...
  openw,lun,'./local/output/GRS_HI_Verification.txt',/GET_LUN
  printf,lun,'CNUM  d_GRS  d_tan  d_ML  d_bar'
  
  hisa_int = fltarr(nind)
  upsilon = FLUX2TAU_BGPS(20.d)
  
  ;; LOOP!  By heaven's sword, LOOP!
  FOR i=0L, nind-1 DO BEGIN
     
     ;; Let the hapless user know where we are in the process...
     IF (i+1) MOD 20 EQ 0 THEN $
        message,'Working object #'+string(i+1,format="(I0)")+' of '+$
                string(nind,format="(I0)"),/inf
     
     ;; Load files for this source
     scnum = string(cnum[i],format="(I04)")
     grs = readfits(post+'grs'+scnum+'.fits',ghdr,/silent)
     hi  = readfits(post+'hisa'+scnum+'.fits',hhdr,/silent)
     lab = readfits(post+'label'+scnum+'.fits',lhdr,/silent)
     smoo = readfits(gpost+scnum+'.fits',shdr,/silent)
     bgps_smoo = readfits(post+'smooth'+scnum+'.fits',bshdr,/silent)
     
     ;; Get the contour plotting coordinates from LABEL image
     plot_wcs_axes,lhdr,10,XC=xc,YC=yc 
     
     ;; Weight Mask -- derived from BGPS Label image
     undefine,wtmask
     wtmask  = aperture_wtmask(lhdr, 40., emaf[i].l, emaf[i].b)
     
     ;; Get HI weighted average inside 40" BGPS beam
     wtmask2     = aperture_wtmask(hhdr, 40., emaf[i].l, emaf[i].b)
     hisa_int[i] = total(wtmask2*hi,/NAN) / total(wtmask2)
     
     tand = !MW.R0*cos(emaf[i].l * !dtor)/1.d3
     
     ;; Check that the postage stamps for GRS & HI exist
     IF total(size(grs,/DIM)) EQ 0 OR total(size(hi,/DIM)) EQ 0 THEN BEGIN
        print,'NO POSTAGE!  ',emaf[i].cnum, emaf[i].l, emaf[i].b, emaf[i].vlsr
     ENDIF ELSE IF plot THEN BEGIN
        
        ;; Open the PostScript file & setup plotting environment
        undefine,grsr,hisar,smoor
        myps,root+scnum+'.eps',xsize=18,ysize=6,ct=0
        
        ;; First, the colorbars
        grsr  = set_plot_range(grs)
        hisar = [-20,20]
        ;;[(set_plot_range(hi))[0],-(set_plot_range(hi))[0]]
        smoor = set_plot_range(smoo)
        
        cgColorbar,/vert,range=grsr,position=[0.05,0.10,0.06,0.90],/yst,$
                   title='GRS Integrated Intensity [K km/s]',charsize=0.7,$
                   format="(F0.2)"
        
        cgColorbar,/top,range=hisar,position=[0.56,0.91,0.86,0.93],/xst,$
                   title='HISA Integrated Intensity (On-Off) [K km/s]',$
                   charsize=0.7,format="(F0.2)",/right
        
        cgColorbar,/vert,range=smoor,position=[0.95,0.10,0.96,0.90],/yst,$
                   title='GLIMPSE Intensity [MJy sr!u-1!n]',$
                   charsize=0.7,format="(F0.2)",/right
        
        multiplot,[3,1],xgap=0.04,mtitle='BGPS #'+scnum
        
        ;;============================
        ;; GRS Panel
        plotmap,grs,ghdr,ct=0,title='GRS',range=grsr,charsize=1.0,$
                min_dpi=100
        
        ;; BGPS flux density contours
        cgLoadct,40,/silent
        cgContour,bgps_smoo,xc,yc,levels=[.025,.05,0.1,0.25,0.5],$
                  c_colors=[1,2,3,4,5]*50,/overplot,thick=1.5,c_charsize=0.7
        
        ;; BGPS 40" aperture & Bolocat source outline
        cgContour,lab,xc,yc,level=[0.5],color=cgColor('Red'),thick=4,$
                  label=0,/overplot
        cgContour,wtmask,xc,yc,levels=[0.5],thick=5,label=0,$
                  color=cgColor('Dodger Blue'),/over,c_linestyle=2
        
        ;; Center of GRS cloud
        cgPlots,emaf[i].grs_l,emaf[i].grs_b,psym=cgSymcat(46),color='Bisque',$
                symsize=2.0
        cgPlots,emaf[i].grs_l,emaf[i].grs_b,psym=cgSymcat(45),color='Opposite',$
                symsize=2.0
        
        multiplot
        
        ;;============================
        ;; HISA Panel
        plotmap,hi,hhdr,ct=0,title='HISA',range=hisar,charsize=1.0,$
                min_dpi=100,XC=hxc,YC=hyc
        cgContour,hi,hxc,hyc,levels=[0.],c_color='Wheat',thick=1.5,$
                  c_charsize=0.7,/over
        
        ;; BGPS flux density contours
        cgLoadct,40,/silent
        cgContour,bgps_smoo,xc,yc,levels=[.025,.05,0.1,0.25,0.5],$
                  c_colors=[1,2,3,4,5]*50,/overplot,thick=1.5,c_charsize=0.7
        
        ;; BGPS 40" aperture & Bolocat source outline
        cgContour,lab,xc,yc,level=[0.5],color=cgColor('Red'),thick=4,$
                  label=0,/overplot
        cgContour,wtmask,xc,yc,levels=[0.5],thick=5,label=0,$
                  color=cgColor('Dodger Blue'),/over,c_linestyle=2
        
        ;; Center of GRS cloud
        cgPlots,emaf[i].grs_l,emaf[i].grs_b,psym=cgSymcat(46),color='Bisque',$
                symsize=2.0
        cgPlots,emaf[i].grs_l,emaf[i].grs_b,psym=cgSymcat(45),color='Opposite',$
                symsize=2.0
        
        multiplot
        
        ;;============================
        ;; GLIMPSE Panel
        plotmap,smoo,shdr,ct=0,title='Smoothed GLIMPSE',range=smoor,$
                charsize=1.0,min_dpi=100
        
        ;; BGPS flux density contours
        cgLoadct,40,/silent
        cgContour,bgps_smoo,xc,yc,levels=[.025,.05,0.1,0.25,0.5],$
                  c_colors=[1,2,3,4,5]*50,/overplot,thick=1.5,c_charsize=0.7
        
        ;; BGPS 40" aperture & Bolocat source outline
        cgContour,lab,xc,yc,level=[0.5],color=cgColor('Red'),thick=4,$
                  label=0,/overplot
        cgContour,wtmask,xc,yc,levels=[0.5],thick=5,label=0,$
                  color=cgColor('Dodger Blue'),/over,c_linestyle=2
        
        ;; Center of GRS cloud
        cgPlots,emaf[i].grs_l,emaf[i].grs_b,psym=cgSymcat(46),color='Bisque',$
                symsize=2.0
        cgPlots,emaf[i].grs_l,emaf[i].grs_b,psym=cgSymcat(45),color='Opposite',$
                symsize=2.0
        
        
        ;; Textual information to be placed at the bottom of the image
        fdata = 1.d - (emaf[i].c / (1 - exp(-upsilon*emaf[i].s40)))
        
        str = 'd!dGRS!n = '+string(emaf[i].grs_dist,format="(F0.2)")+' kpc   '+$
              'd!dtan!n = '+string(tand,format="(F0.2)")            +' kpc   '+$
              'd!dML!n = ' +string(dml[i,0],format="(F0.2)")        +' kpc   '+$
              'd!dbar!n = '+string(dbar[i,0],format="(F0.2)")       +' kpc   '+$
              'V!dGRS!n = '+string(emaf[i].grs_v,format="(F0.1)")+' km/s   '+$
              'V!dBGPS!n = '+string(emaf[i].vlsr,format="(F0.1)")+' km/s   '+$
              'C = '+string(emaf[i].c,format="(F0.3)")+'   '+$
              'S!d1.1!n = '+string(emaf[i].s40,format="(F0.2)")+' Jy   '+$
              'f!ddata,20!n = '+string(fdata,format="(F0.2)")
        cgText,/norm,0.5,0.05,str,align=0.5,charsize=1.0
        
        myps,/done,/mp

     ENDIF
     
     ;; Table row to be written to text file
     str = scnum+'  '+string(emaf[i].grs_dist,format="(F0.2)")+'  '+$
           string(tand,format="(F0.2)")+'  '+$
           string(dml[i,0],format="(F0.2)")+'  '+$
           string(dbar[i,0],format="(F0.2)")
     printf,lun,str
     
  ENDFOR
  
  close,lun
  free_lun,lun
  
  ;; Make a plot showing HISA vs GRS distance
  myps,'./emaf_paper/plots/hisa_vs_grsdist.eps',xsize=10,ysize=4.5
  multiplot,[2,1],xgap=0.04
  
  plot,emaf.grs_dist,hisa_int,psym=8,xtit='d!dGRS!n [kpc]',$
       ytit='I!dHISA!n [K km/s]'
  
  linfit_plot,emaf.grs_dist,hisa_int,linestyle=2,color=cgColor('Red'),thick=2,$
              /legend,linsize=0.6,charsize=0.8
  
  multiplot,/doxaxis,/doyaxis
  
  tand = !MW.R0*cos(emaf.l * !dtor)/1.d3
  
  plot,emaf.grs_dist/tand,hisa_int,psym=8,xtit='d!dGRS!n / d!dtan!n',$
       ytit='I!dHISA!n [K km/s]'
  
  linfit_plot,emaf.grs_dist/tand,hisa_int,linestyle=2,color=cgColor('Red'),$
              thick=2,/legend,linsize=0.6,charsize=0.8
  
  plots,1.1,!y.crange[0]
  plots,1.1,-5,/continue,linestyle=3,color=cgColor('YGB5'),thick=2
  plots,!x.crange[1],-5,/continue,linestyle=3,color=cgColor('YGB5'),thick=2
  
  myps,/done,/mp
  
  
  wtfind = WHERE(emaf.grs_dist/tand GE 1.1 AND hisa_int LE -5, nwtf)
  print,nwtf
  IF nwtf NE 0 THEN print,emaf[wtfind].cnum
  
END
